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I In this paper, three-dimensional (3D) multi-relaxation time (MRT) lattice-Boltzmann 

I (LB) models for multiphase flow are presented. In contrast to the Bhatnagar-Gross- 

' Krook (BGK) model, a widely employed kinetic model, in MRT models the rates 

T— I ■ of relaxation processes owing to collisions of particle populations may be indepen- 

\0 ; 

■ dently adjusted. As a result, the MRT models offer a significant improvement in nu- 
c/3 \ 

. merical stability of the LB method for simulating fluids with lower viscosities. We 
*c/3 ■ 

^~>. show through the Chapman-Enskog multiscale analysis that the continuum limit 

behavior of 3D MRT LB models corresponds to that of the macroscopic dynami- 

^ \ cal equations for multiphase flow. We extend the 3D MRT LB models developed 

H ' 

' to represent multiphase flow with reduced compressibility effects. The multiphase 

models are evaluated by verifying the Laplace- Young relation for static drops and 
the frequency of oscillations of drops. The results show satisfactory agreement with 
available data and significant gains in numerical stability. 
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1 Introduction 



In recent years, computational methods based on the lattice- Boltzmann equa- 
tion (LBE) have attracted much attention. They are based on the paradigm 
of simulating complex emergent physical phenomena by employing minimal 
discrete kinetic models that represent the interactions and spatial and tem- 
poral evolution of quasi-particles on a lattice [1,2]. Originally developed to 
overcome certain drawbacks such as the presence of statistical noise and the 
lack of Gahlean invariance of lattice-gas automaton (LGA) [3], the lattice- 
Boltzmann equation (LBE) [4] has undergone a number of further refinements. 
They include enhanced representation of collisions of particle populations 
through a relaxation process [5] which was further simplified by employing 
the Bhatnagar-Gross-Krook (BGK) approximation [6] later [7,8]. Also, its for- 
mal connection to the Boltzmann equation, developed in the framework of 
non-equilibrium statistical mechanics, was established [9]. 

The LBE has the beneficial feature that by incorporating physics at scales 
smaller than macroscopic scales, complex fluid flows such as multiphase flows 
can be simulated [10]. In particular, phase segregation and interfacial fluid 
dynamics can be simulated by incorporating inter-particle potentials [11], con- 
cepts based on free energy [12] or the kinetic theory of dense fluids [13,14,15]. 
The inter-particle potential approach in Ref. [11], is an earher LBE approach 
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for multiphase flows and is based on non-local pseudo-potentials between par- 
ticle populations. On the other hand, free-energy based methods are derived 
from thermodynamic principles that naturally provide interfacial multiphase 
flow physics [12]. The LBE approaches based on kinetic theory of dense fluids 
represent physics based on mean-fleld interaction forces and exclusion vol- 
ume effects [13,14,15]. More recently, LBE models have also been developed 
to handle high density ratio problems, which also overcome some of the other 
limitations of earlier approaches [16,17]. In general, the approaches above em- 
ployed the BGK approximation to represent the collision process in the LBE 
and some of their applications in 3D are presented in Ref. [18]. 

It is well known that the BGK model, a single-relaxation time model, of- 
ten results in numerical instability when fluids with relatively low viscosities 
are simulated [19]. The instability problems may be compounded in three- 
dimensional (3D) flows when physics may not be adequately resolved owing 
to computational constraints. To address this limitation with the standard LB 
models, several approaches have been proposed. In one approach, known as 
the entropic lattice Boltzmann method(ELBE), the equilibrium distribution 
function is deflned in such a way that it minimizes a certain convex function, 
a Lyapunov functional designated as the i7— function, under the constraint 
of local conservation laws [20,21]. As a result, it ensures the positivity of the 
distribution function of particle populations and thus improves the numerical 
stability. While this approach is endowed with elegant and desirable physical 
features, its numerical accuracy is not established and it would incur relatively 
heavy computational overhead [22]. 

Alternatives to the BGK model have been proposed to improve numerical 
stability. In the multi- relaxation time (MRT) method [23], by choosing dif- 
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ferent and carefully separated time scales to represent changes in the various 
physical processes due to collisions, the stability of the LBE can be signifi- 
cantly improved [19]. Another interesting approach is the fractional time step 
method [25]. In this approach, stability is improved by considering fractional 
propagation of particle populations by reducing time step, which increases 
computational time. While this approach was employed for the BGK model, 
the underlying idea is not limited to a particular collision model. Indeed such 
an approach could be employed in conjunction with the MRT model to fur- 
ther improve stability. Yet another possibility for improving stability is based 
on the so-called two-relaxation time (TRT) models (e.g. Ref. [26]). These are 
variants of the MRT models and can be considered as hmiting cases of the 
more general MRT models. 

In addition to the computational advantage of significantly improved stabil- 
ity, the MRT models also have physical advantages in that they are flexible 
enough to incorporate additional physics that cannot be naturally represented 
by the models based on the BGK approximation. In contrast to the BGK mod- 
els, MRT models deal with the moments of the distribution functions, such as 
momentum and viscous stresses directly. This moment representation provides 
a natural and convenient way to express various relaxation processes due to 
collisions, which often occur at different time scales. Since the BGK model 
can represent only a single relaxation time for all processes, it cannot nat- 
urally represent physics in certain complex fluids. For example, as discussed 
in Ref. [27], BGK models cannot represent all the essential physics in vis- 
coelastic flows in 3D. Moreover, MRT models make it possible to incorporate 
appropriate acoustic and thermal properties through adjustable Prandtl num- 
bers for simulating thermo- hydrodynamics by employing a hybrid LBE/finite- 
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difference model [28]. BGK-based models are not able to do so naturally. Also, 
in a more recent work, MRT models have been shown to be capable of han- 
dhng more general forms of diffusional transport than possible with the BGK 
model [30]. In addition to the capability of dealing with genuinely anisotropic 
diffusion problems, their MRT formulation is able to suppress directional ar- 
tifacts arising due to discrete lattice effects even in isotropic problems. Such 
features are possible with MRT models because they lend themselves readily to 
constructing and controlling various models to evolve at rates consistent with 
the dynamics of physical phenomena of interest. Indeed, approaches based on 
the MRT representation are inspired by the moment method developed in the 
seminal work of Maxwell [31] which was further developed by Grad [32]. 

While prior works have demonstrated the computational advantages of the 
MRT models for single-phase flows in 2D by Lallemand and Luo [19] and in 
3D by d'Humieres et al. [24] and for turbulence modeling [29], it has not been 
shown for multiphase flows in 3D in previous works. Multiphase flows involve 
additional physical complexity as a result of interfacial physics involved - i.e., 
phase segregation and surface tension effects. In this case, the accuracy of the 
numerical discretization of the source terms representing interfacial physics 
becomes an important consideration and needs to be incorporated carefully 
into the framework of the MRT model. These terms should be modeled in 
a way that, when we establish their relation with the various moments, the 
dynamical equations in the asymptotic limit should correspond to the desired 
macroscopic behavior for multiphase flow. These elements were not considered 
in single-phase MRT models. We introduce a second-order discretization of 
these source terms and avoid implicitness through a transformation which, to 
our best knowledge, is applied to MRT models for multiphase flows in 3D for 
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the first time in this work. Recently, a MRT model for 2D multiphase problems 
was developed [34]. There are significant differences in the development and 
implementation of 2D and 3D MRT models. These will be discussed below. 
Furthermore, realistic multiphase flows are 3D in nature and in this work, 
we develop a 3D MRT model for such problems by employing a model for 
interfacial physics based on the kinetic theory of dense fluids [13,14,15]. 

The speciflc contributions of this work will now be stated. We develop a 
3D MRT model incorporating phase segregation and surface tension effects 
through source terms discretized by using trapezoidal rule integration and 
simplifled through a transformation to achieve an effectively exphcit 3D MRT 
model for multiphase flows. Since the underlying lattice structure for 2D and 
3D models are different, the moment basis for the corresponding MRT mod- 
els are different, and is more complicated in 3D. We present the theoretical 
developments based on the 3D moment basis for multiphase flow which, to 
our best knowledge, has not been considered in prior work. We derive the 
continuum equations for multiphase flow from the 3D MRT model through a 
Chapman- Enskog analysis [33], provide the dynamical relations between the 
various moments and forcing terms representing interfacial physics including 
surface tension forces, develop the relationships between the gradients of mo- 
mentum flelds in terms of the non-equilibrium parts of certain moments and 
their corresponding relaxation times, and explicitly derive the relationship be- 
tween the transport coefficients for the fluid flow and appropriate relaxation 
times. To improve the stability of the approach, the model is then transformed 
in such a way that compressibility effects are reduced. Furthermore, we evalu- 
ate the accuracy and gains in stability of the MRT model for some canonical 
multiphase problems in 3D. 
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As the MRT models are endowed with greater computational stability and 
potential to incorporate additional physics, there is considerable interest in 
their apphcations to multiphase problems: The 2D MRT model developed 
by McCracken and Abraham [34] extended for axisymmetric problems using 
the axisymmetric LB model [35] has been employed to study the physics of 
break up of liquid jets [36]. This axisymmetric MRT model and the 3D MRT 
model developed in this work have also been employed to study the physics 
of head-on as well as off-center binary drop coUisions, respectively [37]. The 
physically inspired MRT approach developed in this paper also provides a 
natural framework for incorporating additional physics such as viscoelastic or 
thermal effects in multiphase flows through an LB model, which are subjects 
of future work. The rest of the paper is organized as follows. In Section 2, the 
3D MRT LBE multiphase models are developed. In Section 3, the macroscopic 
dynamical equations of these models are derived by using a Chapman- Enskog 
multiscale analysis. Section 4 transforms the model to simulate multiphase 
flow with reduced compressibility effects. In Section 5, the model is applied to 
benchmark problems to evaluate its accuracy and gains in stability. Finally, 
the paper closes with summary and conclusions in Section 6. 



2 3D MRT LBE Model for Multiphase Flow 



We develop a MRT model for multiphase flows in which the underlying inter- 
facial physics is based on the kinetic theory of dense fluids. In particular, we 
consider that the particle populations representing the dense fluids experience 
mean-field interaction forces and respect Enskog effects for dense non-ideal 
fluids [13,14]. The effect of collisions of particle populations is represented 



7 



though a generahzed relaxation process in which the distribution functions 
for discrete velocity directions approach their corresponding local equilibrium 
values at characteristic time scales given in terms of a generalized collision or 
scattering matrix. 

In the following, we consider subscripts with Greek symbols for particle ve- 
locity directions and Latin symbols for Cartesian components of spatial direc- 
tions. We assume summation convention for repeated indices for the compo- 
nents of spatial directions. In addition, unless otherwise stated, we follow the 
convention that vectors corresponding to three-dimensional position space are 
represented by non-capitalized symbols with arrowheads; the vectors with a 
particle velocity basis are denoted by non-capitalized boldface symbols; the 
square matrices constructed from the particle velocity basis are represented by 
non-boldface capitalized symbols. We consider the following MRT LBE with 
a source term that gives rise to phase segregation and surface tension effects: 

/a(^ + et5t,t + St) - fa{'^,t) = ^a\ix,t) + 



Here, fa is the discrete single-particle distribution function, corresponding to 
the particle velocity, e^, where a is the velocity direction. The Cartesian com- 
ponent of the particle velocity c, is given by c = S^/St, where is the lattice 
spacing and St is the time step. The left hand side (LHS) of this equation rep- 
resents the change in the distribution function as particle populations advect 
from one lattice node to its adjacent one along the characteristic direction 
represented by the discrete lattice velocity direction (see Figs. 1 and 2). 
On the other hand, the right hand side (RHS) represents the effect of particle 
collisions and force interactions. 





8 



The form of the coUision term that incorporates MRT coUision processes 
will be discussed below. The source term Sa which models the interfacial 
physics in multiphase flow needs to be accurately represented. In Eq. (1), we 
considered a second-order trapezoidal rule discretization of this term. It may 
be written as [13] 

Sa = (2) 

where Fj is the interaction force term that models the interfacial physics, 
while Fgxt,j corresponds to external or imposed forces such as gravity. F^ is 
modeled as a function of density following the work of van der Waals [38]. 
The exclusion- volume effect of Enskog [33] is also incorporated to account for 
increase in collision probability due to the increase in the density of non-ideal 
fluids. These features account for the phase segregation and surface tension 
effects. For more details, the reader is referred to Ref. [13]. In Eq. (2), /^^'*^ is 
the local discrete Maxwellian and its functional expression will be presented 
below. Effectively, Fj may be written as 

Fi^-d^^l; + F,j. (3) 

In Eq. (3), Fgj represents the surface tension force and is related to the density 
p and its gradients by 

F,,, = KpdjV'p, (4) 

where k is a surface tension parameter. It is related to the surface tension a 
of the fluid by the equation [40] 
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where n is the direction normal to the interface. Thus, the surface tension is a 
function of both the parameter k and the density profile across the interface. 

The term ip in Eq. (3) refers to the non-ideal part of the equation of state 



is employed where 7 = bp /A. The parameter a is related to the inter-particle 
pair-wise potential and b is related to the effective diameter d of the particle, 
and the mass m of a single particle, by 6 = 2TTd^/3m. ip assumes an important 
role in determining phase segregation [13]. The Carnahan-Starling-van der 
Waals EOS has aP — 1/p — T curve, in which dP/dp < for certain range 
of values of p, when the fluid temperature is below its critical value. This 
part of the curve represents an unstable physical situation and is the driving 
mechanism responsible for keeping the phases segregated and in maintaining 
a self-generated interface that is diffuse with a thickness of about 3-4 lattice 
grid points. 

The local discrete MaxweUian /^^'^(p, it) in Eq. (2) is obtained from a trun- 
cated expansion in terms of fluid velocity it of its continuous version, after 
approximating it for discrete particle velocities [9]. As a result, it is a func- 
tion of the fluid densities and velocities and is given by 



(EOS) 



ll)[p)^P- pRT. 



(6) 



In this work, the Carnahan-Starling-van der Waals EOS [39] , 




(7) 
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i. 



Fig. 1. D3Q15 lattice. 



where Wa is the weighting factor. In this representation of the local Maxwellian, 
the factor RT is related to the speed of sound Cg of the model through RT — 
c^, where Cg — 1/V^c. For the three-dimensional, fifteen-velocity (D3Q15) 
model [7], shown in Fig. 1, the weighting factors become 



a — 1 



h a = 2, •••,7 



(9) 



and for the three-dimensional, nineteen- velocity (D3Q19) model [7], shown in 
Fig. 2, we have 
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Fig. 2. D3Q19 lattice. 



18 ^ — 2, • • • , 7 



2g 8 ^ * * * ^ 1 9 • 



The corresponding particle velocity directions are 



(0,0,0) 



a — 1 



:±l,0,0),(0,±l,0),(0,0,±l)a = 2. 



(d=l,±l,±l) 



a = 8, 
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and 



(0,0,0) 



a = 1 



(±1,0,0),(0,±1,0),(0,0,±1) a = 2, ■■■,7 



:±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1) a = 8, • • • , 19, 



(12) 



respectively. 

We express the effect of collisions in Eq. (1) as a relaxation process through 
a multi-relaxation time (MRT) model for the collision term as shown be- 
low [19,23,24]: 



/3 



(13) 



where A^p is the component of the collision or the scattering matrix A. Here, 
the equilibrium distribution is an appropriate function of the conserved 
moments such as the density and momentum and is, in general, not necessarily 
the local discrete Maxwellian given by Eq. (8). The hydrodynamic field vari- 
ables are obtained by taking the kinetic moments of the distribution functions 
as 



a 



(14) 

(15) 



Equation (1) is implicit. For computational convenience, it can be made ex- 
plicit if we introduce the transformation [13] 



fa — fa 2*^"^*" 



(16) 
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As a result, Eq. (1) is replaced by the following effectively explicit scheme: 

faC^ + e^St, t + 6t)- faC^: t) = - ^ A^/J (//j - /|'') | + 

(jocfi - 2^«/3) Sp\{x,t)St, (17) 
where lap is the component of the identity matrix I. 

Next, we construct appropriate sets of linearly independent moments from 
the distribution functions in velocity space. Since moments of the distribution 
function, such as the momentum and viscous stresses, directly represent phys- 
ical quantities, the moment representation offers a natural and convenient way 
to express the relaxation processes due to collisions. One particular advantage 
of this representation is that the time scales of the various processes repre- 
sented in terms of moments can be controlled independently. The moments are 
constructed from the distribution function through a transformation matrix 
T comprising a linearly independent set of vectors, i.e. 

/ = Tf, (18) 



where 



/ — fi, /2, fs, fi, fb, /e, /?, fs, /g, fio, fll, fl2, fl3, /l4, fib 



(19) 



is the vector representing the particle distribution functions for the D3Q15 
model and 



/ — /l; /2) /si /4j /si /ei /ti /s; fd: flO, fll: fl2, flS: 

T 



fl4, fib, fw, fn, flS, fl9 



(20) 



for the D3Q19 model. 



14 



In Eqs. (18) - (20), the superscript 'T' represents transpose of a matrix and 
the vector / represents the moments 

S — P-i^-i^ J jxj Qxj jyj Qyj jzj Qzj ^Pxxj Pwwj Pxyj Pyzj Pzxj TTT'xyz i (^1) 



and 



f — 1 jx-i Qxi jyi Qyi jz-i Qzi ^Pxxi ^'^xxi Pwwi 



'^■wwiPxyiPyz-iPzx-i'^x-i'^yi'^z\ i (^^) 

for the D3Q15 and D3Q19 models, respectively. For the D3Q15 model, e and 
represent kinetic energy that is independent of density and square of energy, 
respectively, j^i jy and are the components of the momentum or mass 
flux {jx — pux,jy — puy,jz — puz), Qx, Qy, Qz ^re the components of the 
energy flux, and p^x, Pxy, Pyz and p^x are the components of the symmetric 
traceless viscous stress tensor. The other two normal components of the viscous 
stress tensor, pyy and Pzz, can be constructed from pxx and p^y,, where p^y, — 
Pyy ~Pzz- The quantity rUxyz is an antisymmetric third-order moment. For the 
D3Q19 model, instead of the moment nixyz-i we have five additional moments: 
37ra;a;, 37r^^, ma;, my and m^. The first two of these moments have the same 
symmetry as the diagonal part of the traceless viscous tensor pjj, while the 
last three vectors are parts of a third rank tensor, with the symmetry of 

ikPmn [24]. 

The underlying principle for the construction of the transformation matrix is 
based on the observation that the collision matrix reduces to a diagonal form 
in an appropriate orthonormal basis, which is obtained by combinations of 
monomials of Cartesian components of the particle velocity directions [23]. 
Following d'Humieres et al. [24], for the D3Q15 model we introduce a trans- 
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formation matrix T, which represents components of the 15 orthogonal basis 
column vectors Vj} that form the moment basis, i.e. 

T = [vi, V2, V3, V4, V5, va, vr, Vs, Vg, Vio, Vu, v^, ^13, Vu, Vi^f , (23) 

where the components of each column vector may be written as 

Via = |e;^|°; V2a = |e^p-2; = KlSle^l"^ - 55|e^p + 32) ; = e^x', 

Vba — 2(^1^"^!^ 13) Cqj;; Vqq. — ^ay'i — 2 13) Cqj/; Vga — ^az'i 

V<d(x — \ (5|e^P — 13) Cq^; t'loa — ^^ax ~ I^Pj Vila — ^ay ~ ^az'i Vl2a — ^ax^ay] 
Vl3a ^ay^az^ ^14a ^ax^azj ^15q; ^ax^ax^az- 

For the D3Q19 model, we have 

T= [Vi, V2, V3, V4, V5, Vq, Vt, Vs, Vg, Viq, Vn, Vi2, Vi3, 

Vi4,*^i5,'yi6, Vi7,'yi8,'yi9]'^ (24) 

based on 19 orthogonal basis column vectors. The components of each of these 
vectors are 

Via = |e^|°; V2a = 19|e^P - 30; v^a = I (21|e^|^ - 53|e^p + 24) ; V4a = eax] 

V5a (5|Ca| 9) Caxi ^6a *^ayj ^7a (5|Ca| 9) Cayj ^8a ^azi Vga 

(5|e^P - 9) Caz] vioa = 3el^ - le^p; vua = (3|e2p - 5) (3e^^ - |e^p) ; vua = 

^ay ^azi Vl3a — (3|Ca| 5) (^G-^y ^az^ i Vl4a — ^ax^ayi ^15a — ^ay^azi ^16a — 
^ax^az) — (^oy ~ ^az) ^axj VigQ. — (e^^ — 6^^) C^y, Vig^ — ((^ax ~ ^ay) ^oiz- 

These matrices are formed by a set of linearly independent orthogonal basis 
vectors (i.e. Va • vp — la^a/s, where la is a normalizing constant) which are 
constructed by a Gram-Schmidt procedure such that they diagonalize the 
collision matrix A i.e. A = TAT~^ [24], where the matrix A in moment 
space is a diagonal matrix. This orthogonal set is built in increasing order of 
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moments and then arranged in increasing order of complexity of the tensorial 
representation of the moments. An example of such a construction with some 
details for a 2D MRT model is given by Bouzidi et al. [41]. As an example, 
for unit lattice spacing and time steps, i.e. c = 1, the simplified form of the 
transformation matrices is given in the Appendix. 

The equilibrium distribution functions / in moment space are related to 
those in velocity space, /^^, through the transformation matrix 



n p<^1 p2,eg ■ eq ■ eq ■ eq o eg eq eq eq eq eq 1^ /r,c-\ 
Hi^ 1^ )Jx)Hx iJyi'iy tJztHz ■> ^yxx-> ywwi yxyi yyzi yzxi '"'xyz l^^/ 



for the D3Q15 model and 



^ p^q p'2,eq ■ eq ■ eq ■ eq n eq o eq eq eq 
1^ y Jxy ^ Jy^ ^ J ZT Hz T '-'fxxi '-"^ wwi I'wwi wwi 

„e9 „e<? „eg eq eq ^eg"''^ 
yxyi Pyzi yzxi '"'x ' '"'y ' '"'z 



(26) 



for the D3Q19 model. Notice that for the conserved or hydrodynamic moments 
corresponding to density and components of momentum, i.e. fp, where /3 = 
1,4,6,8, their equilibrium distributions in moment space are also the same, 
i.e. = fp, where (3 — 1,4,6,8 since the collision process does not alter 
hydrodynamic moments. 

The expression for the equilibrium distribution functions of the non-conserved 
or kinetic moments, which are in turn algebraic functions of the conserved 
moments and obtained by optimizing isotropy and Galilean invariance, are 
given by [24] 

psq _ _„ I i^x+jy+jz) . 2,eq _ ^ _ r (jI+J^+jI) . eq _ _7 ■ . eq _ _7 ■ . eq _ 

— fJ ^ p ,c — fj o p 1 Hx ~ 3^2:1 Hy — r^Jyi Hz ~ 

— In ■ n^l — 1 [^i^-Qy+iz)] . eq _ \^r2±\.. ^eq _ U3y . eq _ M±- „eg _ i^U . 
•iJzi fxx 3 p 1 fww p 1 Pxy p ' Pyz p ' Pxz p ' 
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m 



eq 
xyz 



for the D3Q15 model and 

ge. = -lip + 19 (^^±11^; e^,eq ^ 3^ _ n M±|+i|l. ^e, ^ _2^-^. ^e, ^ _|^-^. 

„eq _ _2„- . „eq' _ 1 . eg _ _l„eg. „eg _ [j^"-?^] . eg _ _l„eg . 

"iJ^i i^xx 3 p 1 '^xx l"xxi cww p 1 "ww 2^ww) 

Pll = '-f'; Pll = Pll = "f; ^7 = O; = O; = 
for the D3Q19 model 

The coUision matrix A in the moment space is given by 

A = diag [si, S3, S4, S5, Sq, s^, Ss, Sg, Sio, Sii, 512, S13, S14, 515]'^ (27) 

for the D3Q15 model and 

A= diag [Si, S2, S3, S4, S5, Sg, Sy, Sg, Sg, SiQ, Sii, S12, Sl3, 

Sl45 Si5, S16, Si7, S18, (28) 

for the D3Q19 model, where the parameters sp represent the inverse of the 
relaxation times of the various moments / in reaching their equilibrium values 

r. 

The variables si, S4, sq and ss are the relaxation parameters corresponding 
to the collision invariants p, j^, jy and j^, respectively. Since the collision 
process conserves these particular moments, the choice of their corresponding 
relaxation times is immaterial. This is because they are each specified to re- 
lax during collisions to their local equilibrium which are actually defined to be 
the corresponding pre-coUision value of the respective quantities in the equilib- 
rium distribution /"^^ (see Eqs. (25) and (26)). So, whatever be the relaxation 
parameters for these quantities, the collision process does not change their 
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values. The relaxation parameters for these collision invariants could be set 
to zero when there is no forcing term Sa in the LBE. However, care needs to 
be exercised in selecting the relaxation parameters of MRT LBE with forcing 
term. This has been pointed out in Ref. [42]. This is because the collision ma- 
trix also influences the forcing term in the effectively explicit MRT LBE (see 
Eq. (17)). To correctly obtain the continuum asymptotic limit of this LBE, i.e. 
the macroscopic hydrodynamical equations for multiphase flow, derived in the 
next section, the relaxation times for the conserved moments should be set to 
non-zero values. For simplicity, we have chosen a value of unity as the relax- 
ation parameter for these moments (si = S4 = sq = ss = 1) that preserves the 
influence of the second-order discretization of the source terms in Eq. (1) (see 
also Ref. [34]). Furthermore, for the D3Q15 model, it will be shown through 
the Chapman-Enskog analysis in Section 3 that the parameter S2 is related 
to the bulk viscosity and Sio through S14 are all related to the kinematic vis- 
cosity. This leaves us with S3, S5, S7, sg, and S15 as free parameters in this 
model. On the other hand, for the D3Q19 model by analogy, the parameter 
S2 is related to the bulk viscosity and sio, S12, Si4 through sie are each related 
to the kinematic viscosity. Then, S3, S5, sy, Sg, Su, Sis, S17 through s^q are the 
free parameters. 

To obtain the LBE model in moment space, the source terms are also trans- 
formed as follows: 

? = Tq, (29) 



where q is the column vector corresponding to the components of the source 
term Sa', ? corresponds to the components of Sa- Now, prc-multiplying Eq. 
(17) by the transformation matrix T, we obtain the 3D MRT model 
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- (//? ~ fp) I (^,*) + Z] ( -^«/3 - ^Aa^ ) -^/Sl (a;,t)(^t- (30) 

The hydrodynamic field variables can be obtained from the distribution func- 
tions in velocity space as 

P=E/- (31) 

a 

pui = E heai + \{Fl + Fe,t,i)St, (32) 

or, more obviously, directly from the components of the moment space vector, 
/, i.e., fi and where (3 — A, 6, 8 for density and momentum, respectively. 



3 Macroscopic Dynamical Equations for Multiphase Flow 



In this section, the macroscopic dynamical equations for the 3D MRT LBE 
multiphase flow model developed above will be derived by employing the 
Chapman- Enskog multiscale analysis [33] for the D3Q15 velocity model. The 
macroscopic equations for other 3D MRT models, such as the D3Q19 model 

discussed in the previous section, can be found in an analogous way. Introduc- 
ing the expansions [43] 



oo 



n=Q 



A„ = 9t„ + Cakdk, 

oo 

n=0 

oo 

n=0 



(33) 
(34) 

(35) 

(36) 
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where e — 5t in Eq. (1), the following equations are obtained as consecutive 
orders of the parameter e: 



.(2) 



(37) 
(38) 

(39) 



In equivalent moment space, obtained by multiplying Eqs.(37)-(39) by the 
transformation matrix T, they become 



0{e') : 9, + Ao (lap - lAa,) = - E ^o^pfl 



(40) 
(41) 

(42) 



where £i — TtaiT ^. 



First let us simphfy the source term Sa in Eq.(38) as 



~n {^az "^z) ~l 7 (^a ' ) ^az 



Fx + 

Fy + 



(43) 



by neglecting terms of the order of 0{Ma?) or higher in its definition, i.e. 
Eq.(2) and Eq. (8). In Eq. (43), F^;, Fy and are the Cartesian components 
of the net force experienced by the particle populations including those that 
lead to phase segregation and surface tension effects and imposed forces (e.g., 

Fx (^x'4^ ~l~ -^5,0; ~l~ Fg^xt,x)- 



Upon substituting Eq. (43), after pre- multiplying by T, in Eq. (41), we get 
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the components of the first-order equations in moment space, 



dtoP + dxjx + dyjy + dj^ = 0, 



dto y-p + • Tj - ^ (djx + dyjy + dj^ 
-S2e^^^ + 2 {F^Ux + FyUy + F^Uz) , 



dto (p-^y - ^ {dxjx + dyjy + dj,) = 
-sae^^^^ - 10 {FxUx + FyUy + F,u,) , 



dtjx + dx r-p + ^ J + 5. f^j + d. 



9^o{-lj.)+dx(-lp + j-^ 



-7 fx + 5^ + 5jI 



+ 



^5Hx 3 



P 



P 



dt. ( -^j.) + dx (m^) a, (-Ip + 1 [5,- - 7,- + 5,- 



JyJz 



-srqy 



w -If 

y 3^2/, 



3' P 



d 



to 
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-Ssiopi'i + 2 (2F,H, - FyUy - F,u,) , (53) 



•2 _ -2 
Jy Jz 



2 

+ 3 [dyJy - dzJz] = 



-Ssiipi'l, + 2 {2FyUy - F,u,) , (54) 



dto [^^^h^ + ^ [^^•^'2/ + ^f-^'^l = -3si2P^^J + F^Uy + FyU^, 
dto {~^y^^ + \ \Pyj^ + ^^-^2/] ^ -3si3p£^ + F^m^ + F^Uy, 



(55) 

(56) 
(57) 



and 



The evolution equations express variations at the to time scale level, the dy- 
namical relationships between various moments constructed from the velocity 
space, their relaxation parameters and the net forces that include those that 
lead to phase segregation and surface tension effects and any external forces 
acting on the particle populations. 

Similarly, the components of the second-order equations in moment space, i.e. 
Eq. (42), can be obtained. Our interest is in the dynamical equations for the 
conserved moments. For brevity, here we express below only the second order 
equations of the conserved moments representing variations at the ti time 
scale. 

dt.p^O, (59) 
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dtjx + ( 3 



1 So 

2 



1 - 2^10^ 



I + 



ft, 



1 - 2^12^ 



1 - 2^14 



0, 



(60) 



5*1 + 



1 Si2 

2 



1 - -S2 

2 



=(1) 



1 
2 



+ 



1 - 2^10^ 



:Si3 



= 0, 



1 - 2^11 



(1) 



+ 



(61) 



and 



dtiiz + dx 



■Si4 



1 - -S2 

2 



e(i) _ i 
2 



1 - 2^10^ 



i Si3 

2 . 
_ 1 

i^a;a; 2 



1 Sii 

2 



trww 



0. (62) 



Now, combining the first- and second-order equations for the conserved mo- 
ments, i.e. Eqs. (44) and (59), Eqs. (47) and (60), Eqs. (49) and (61) and Eqs. 

(51) and (62), by using ft = ft^ + eftj, we get 



ftp + dxjx + dyjy + d^j^ = 0, 



(63) 



OtJx + dx\-P+ -Jx + 



1 S2 

2 



1 - 2^10^ 



Px 



Jxjy 



+ e 



1 - 2^12^ 



dz h e 



1 - -.14 



(64) 



ftj. + 5.(^^H-6 



1 Si2 

2 



1-2^2 



e(^) - e- 



1 - 2^10^ 
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1 - 



+ h e 



1 - 2^13 



(65) 



and 



dtjz + d, 

JyJz 



Jxjz , 

X 1 \- e 



1 6" 14 

2 



Pzx 



+ 



1 



+ e 



P 



1 - 



,(1) 



1 

-e 
2 



1 - -.10^ 



Irxx 



1 

-e 
2 



1 fill 

2 



Irww 



Fz, (66) 



respectively. In these equations, e^^^ and p^^j, p^^;) PxJ' ^^a;'' 
knowns to be determined. Writing expressions for these terms from Eqs. (45) 
and Eqs. (53)- (57), respectively, and employing the continuity and momen- 
tum equations, i.e., Eqs. (44), (47), (49), (51), and neglecting terms of order 
0{Ma^) or higher, we get 



2 2 1 ^ 

3^^^ = -o — i^xix + dyjy + dzjz) = — " / 

2 1 

- — {2d^j^ - dyjy - d^j^) 

y -sio 

6 Sii 
1 1 



e — e 



eq 



fxx 



pW : 

fww 



Pxx Pxx 1 



z) ^ fww 



WWl 



fxy 



t^yz 



pH) 
yzx 



-J -512 

— {dyjz + djy) ^ Py^ - PI% 

- o — i^xjz + d,j^) fti p^^ - pII- 

6 Su 



(67) 
(68) 
(69) 
(70) 
(71) 
(72) 



These equations represent the various gradients of the mass flux fields or 
momentum in terms of non-equilibrium parts of certain moments and their 
corresponding relaxation times. 



Using Eqs. (67)-(72) in Eqs. (64)-(66), the momentum equations simplify to 
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dtjx + d^c 



+ dy 



-9x ( ) + + <9a; 



JyJx 
. P 

2u 



Jzjx \ 

. p )~ 



J + 



dy {l^ {dxjy + dyjx)) + Qz {dxjz + d^jx)) , 



(73) 



dtjy + 



Jxjy 



a, 



Jzjy 



-dy +Fy + d^ (i/ {djy + dyj^)) + 



dy ( 2v 



dyjy -y + ("^ ■ 7^ + dz {dyjz + dz3y)) , 



and 



dtjz + dx 



Jxjz 



d,. 



JyJz 



-dz [:^p] + Fz + d^{iy {d^jz + dzjx)) + 



9y (i^ {dyjz + dzjy)) + dz(2u 



dzjz - • 7 



(74) 



(75) 



where C and i/ are the kinematic bulk and shear viscosities repectively. These 
are related to the relaxation parameters of the energy and the stress tensor 
moments by 



^=^^-^)5t,/3 = 10,ll,12,13,14. 



(76) 
(77) 



Notice that some of the relaxation parameters {siq — su — ■ ■ • — su) should 
be equal to one another to maintain the isotropy of the viscous stress tensor. 
From Eqs. (76) and (77), we see that the bulk and shear viscosities can be 
chosen independently. In particular, this allows us to choose a bulk viscosity 
larger than the shear viscosity, so that acoustic modes are attenuated quickly. 
This is helpful in certain physical situations and could also aid in improving 
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numerical stability. Finally, substituting for the forcing terms, i.e. Fj = —diip+ 
Fg^i + Fext,i-i and recognizing that the pressure is given by P = -0 + l/3p, the 
macroscopic dynamical equations for the conserved moments, i.e. density and 
momentum, are given by 



p ^dtUx + It ■ "^Ux 

P \dtUy + it ■ "^Uy 

p ^dfU;, + it ■ "^u^ 



= 0, (78) 

= d,P + F,,, + Fe,,,, + d^al, + dyaly + d^al,, (79) 

-- dyP + Fs,y + Fe,*,, + d^al, + dyal^ + (80) 

-- d,P + F,,, + Fe,t,, + + dyaly + (81) 



where a^, is the deviatoric part of the viscous stress tensor. 



{dj [pui] + di (puj)) - ■ (pit) 6ij + • (pit) Sij, (82) 



F5 j is the surface tension force given in Eq. (4), Fg^t^i is the external force, and 
i,j e {x, y, z}. The equations (78)-(81), indeed, correspond to the macroscopic 
equations for multiphase flow derived from kinetic theory [45] . 



4 3D MRT LBE Model for Multiphase Flow with Reduced Com- 
pressibility Effects 

Multiphase LBE models based on inter-particle interactions [13] face diffi- 
culties for fluids far from the critical point and/or in the presence of external 
forces [44] . This difficulty is related to the calculation of the inter-particle force 
involving the term djip, in the model which becomes quite large across inter- 
faces. As a remedy. He and co-workers [14] developed a suitable transformation 
of the distribution function f^, by invoking the incompressibility condition of 
the fluid, to Qa that determines the hydrodynamic fields, i.e. pressure and ve- 
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locity fields. A separate distribution function is also constructed to capture 
the interface through an order parameter, which is different from density. The 
introduction of these features helps in reducing compressibility effects in mul- 
tiphase problems. In this section, we apply this idea to the 3D MRT model 
discussed in Section 2. Following He et al. [14], we replace the distribution 
function by another distribution function through the transformation 

g^ = URT + i^ip Y^ K (83) 

By considering the fluid to be incompressible, i.e. 

^V'(p) = (9t + iifeafe)V'(p) = 0, (84) 

and using Eqs. (83) and (16), Eq. (1) becomes 
9a{'^ + e^Su t + St)- ^a(^, t) = 

- J2 ^o^P fe ~ df) I {x,t) + X] ( 4/3 - -^^ap ) Sgp I (85) 

where 

fj = rjRT + i^ipY-^^^^. (86) 



The corresponding source terms are 



(87) 



VP P 

Equations (85)- (87) are numerically more stable than the original formulation 
since the term djip is multiplied by a factor proportional to the Mach number 
0{Ma), instead of being 0(1) when the "incompressible" transformation is 
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not employed. As a result, this term becomes smaller in the incompressible 
limit. Now, applying the transformation matrix to the above system, we get 



In this new framework, we still need to introduce an order parameter to cap- 
ture interfaces. Here, we employ a function referred to henceforth as the 
index function, in place of the density as the order parameter. The evolution 
equation of the distribution function, whose emergent dynamics govern the 
index function, must have the term responsible for maintaining phase seg- 
regation and mass conservation. In this regard, we use Eqs. (1) and (2), by 
keeping the term involving dj%l) and dropping the rest as they play no role in 
mass conservation. In addition, the density is replaced by the index function 
in these equations. Hence, the evolution of the distribution function for the 
index function is given by 




(88) 



/a(^ + 



e^5t,t + 5t) - /a(x',i) = -ZA„/3 Up - -/I") \{x,t) + 




(89) 



where 




a 



pRT 



Ja 



feq,M 



{P, ^)- 



(90) 



In moment space, this becomes 



/a(^ + e^St, t + 6t)- /^( x", t) 




(91) 
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It follows that hydrodynamical variables, such as the pressure and fluid veloc- 
ity, can be obtained by taking appropriate kinetic moments of the distribution 
function g^, i.e. 

P-E9a-lu,d,^{p), (92) 
pRTu, = QaSai + ^RT (F,,, + F,,t,,) St. (93) 

The index function is obtained from the distribution function by taking the 
zeroth kinetic moment, i.e. 

= E/- (94) 



The density is obtained from the index function through linear interpolation, 
i.e. 

p(</')=Pl + /^(p^/-Pl), (95) 

where pl and pu are densities of light and heavy fluids, respectively, and 0l 
and refer to the minimum and maximum values of the index function, re- 
spectively. These limits of the index function are determined from Maxwell's 
equal area construction [38] applied to the function ip{4>) + 4>RT. If the viscosi- 
ties are different in different phases, the corresponding relaxation parameters 
(sii through Si4, which are equal to one another) in each phase are obtained 
from Eq. (77) and their variation across interfaces are determined through a 
linear interpolation similar to that for the density as shown in Eq. (95). 

The numerical implementation of this MRT model is as follows: In the collision 
step, the evolution of distribution functions are computed in moment space 
{f and g) through direct relaxation of moments at rates determined by the 
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diagonal collision matrix A with appropriate updates from the source terms. 
After transforming the post-collision moments back to the velocity space, the 
streaming step is performed in velocity space (/ and g). After these two com- 
putational steps, the various macroscopic fields are updated. This allows an 
efficient implementation of the MRT model for multiphase problems. 

In terms of the computational memory and time resources, the MRT model 
requires a negligible increase in memory and a moderate increase in compu- 
tational time, when compared to the BGK model, if appropriate optimization 
strategies are taken. Depending on the lattice velocity model, the transforma- 
tion matrix T is either a 15 x 15 or 19 x 19 matrix. These matrices are the same 
throughout the domain and so need not be stored at every lattice site, thus 
contributing to negligible increase in memory. The corresponding relaxation 
time collision matrices A in moment space which are diagonal matrices with 
15 or 19 elements, respectively, also do not contribute to any significant in- 
crease in memory requirement. The equilibrium distributions in moment space 
depend only on the local macroscopic field variables and can be computed lo- 
cally without requiring additional storage. The memory requirement for the 
distribution functions for the collision and streaming steps in the MRT model 
is about the same as that required for the corresponding BGK version. 

As noted in Ref . [24] , suitable code optimization techniques should be apphed 
to avoid substantial increase in computational time due to the use of the MRT 
model. First, direct matrix computations, which can lead to increased com- 
putational overhead, should not be carried out in the transformation between 
velocity and moment spaces. Instead, the transformation between the com- 
ponents of the distribution functions in velocity space fa and the moments 
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fp should be carried out explicitly using the expressions obtained from the 
transformation matrix, which is an integer matrix with several common ele- 
ments and many zeroes. In these calculations, all the common sub-expressions 
should be computed just once. Second, the computational algorithm should 
be carried out as discussed in the previous paragraph. With such optimization 
strategies, the MRT model only incurs a moderate increase in computational 
time, about 35% more than the BGK model for the cases discussed in the next 
section, but with significantly improved stability. 

5 Results and Discussion 

The 3D MRT model developed in the previous section will now be evaluated 
for some benchmark multiphase problems. Unless otherwise specified, we ex- 
press the results in lattice units, i.e., the velocities are scaled by the particle 
velocity c and the distance by the lattice spacing S^. We first consider a static 
multiphase problem, namely, the verification of the well-known Laplace- Young 
relation for a stationary 3D drop. This relation expresses the force balance be- 
tween the excess pressure inside the drop and the surface tension force. If AP 
is the pressure difference between the inside and outside of a drop of radius 
Rd with surface tension a, we have AP = 2a /Rd- 

We first simulate this problem with a D3Q15 lattice by considering drops of 
three different radii 15, 21 and 27. The corresponding 3D domain is discretized 
by 41 X 41 X 41, 61 X 61 X 61, and 81 x 81 x 81 lattice sites, respectively. Periodic 
boundary conditions are considered in all directions. We simulate drops in a 
gas with an equal shear kinematic viscosity of 1.0 x 10~^ for the gas and 
liquid and a density ratio of 4. We have chosen equal kinematic viscosities for 
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both the phases just for simphcity. The MRT model developed in this work is 
not limited to choosing only identical viscosities for both the phases. Indeed, 
as noted in the previous section, when there is a viscosity mismatch, the 
corresponding relaxation parameters will be different in each phase and their 
values across interfaces can be obtained by linear interpolation of viscosities 
through the index function 0. It is a known limitation that the underlying 
LB interfacial physical model [13,14] is stable for only low density ratios. The 
MRT model has relatively little influence on this limitation. However, as will 
be shown below, it can significantly increase the stability at lower viscosities. 

It may be noted that when the BGK model is employed to simulate multiphase 
flow with the viscosity used above, the computations become unstable. As 
discussed in Section 2, we choose the relaxation parameters Si, S4, Sq and Sg 
to be 1.0. The parameters Sa, where a = 10, 11, 12, 13, 14, as shown in Eq. 
(77), are obtained from the shear kinematic viscosity. On the other hand, the 
bulk kinematic viscosity is related to S2 (see Eq.(76)). We consider S2 = 1-0. 
Thus, the time scales for the shear and bulk viscosities are different. The 
remaining free parameters S3, S5, sy, Sg and S15, whose values have no effect 
on hydrodynamics, are chosen to be equal to 1.0 for simplicity. 

Figure 3 shows the computed pressure difference across the drop interface as a 
function of inverse drop radii for three values of the surface tension parameter: 
K, — 0.08, 0.10 and 0.12. The lines correspond to a linear fit for the computed 
data. We find that for these selected parameters the computed pressure dif- 
ference agrees with the Laplace- Young relation prediction to within 8%. To 
check the vahdity of the MRT lattice model with a larger velocity set, i.e. the 
D3Q19 model, we performed some selected computations, which are shown 
as open symbols in the same figure. The D3Q15 and D3Q19 models yield re- 
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Fig. 3. Pressure difference across a drop as a function of radius for different values 
of the surface tension parameter (k); Bold symbols employ the D3Q15 model, open 
symbol the D3Q19 model. 

suits that have negUgible difference between them. Henceforth, for comparison 
purposes we will only consider results obtained with the D3Q15 lattice. The 
density profile across the drop interface is found to be characterized well with 
no anisotropic effects. As with other LB models, and with methods based on 
the direct solution of Navier-stokes equations for multiphase flows, velocity 
currents around the interfaces are observed. These are generally small and is 
found to be proportional to the surface tension parameter k. 

Next, we consider a dynamical problem, namely, the oscillation of a liquid drop 
immersed in a gas. We employ the analytical solution of Miller and Scriven 
(1968) [46] for comparison with the computed time periods. According to 
Ref . [46] , the frequency of the n*'* mode of oscillation of a drop is given by 



LUn = CO,, 




(96) 
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where a;„ is the angular response frequency, and uu* is Lamb's natural resonance 
frequency expressed as [47] 




n(n+ l)(n- l)(n + 2) 
Rl [npg + {n+ l)pi] 



(97) 



Rd is the equilibrium radius of the drop, a is the interfacial surface tension, 
and pi and Pg are the densities of the two fluids. The parameter x is given by 



where pi and pg are the dynamic viscosities of the two fluids. Here, the sub- 
scripts g and / refer to the gas and liquid phases, respectively. We consider the 
second mode of oscillation, i.e. n — 2, and analytical expression for the time 
period Tanai is obtained from Eq. (96) through Tanai = '2n/u2- 

The initial computational setup consists of a prolate spheroid with a minimum 
and maximum radii of 11 and 15, respectively, and placed in the center of a 
domain discretized by41x41x41 lattices. We consider a density ratio of 4 
for different values of shear viscosity and surface tension parameters. Figure 4 
shows the drop configurations at different times for u — 1.6667 x 10~^ and 
K — 0.10. Att = 100, the drop has a prolate shape which temporarily assumes 
an approximate spherical shape at t = 500. At t = 700, notice that the 
drop becomes an oblate spheroid and at the longer time oi t — 6000, the 
drop assumes its equilibrium spherical configuration. Computations are also 
performed by reducing the shear viscosity by 2.5 and 5 as a function of time 
from the above value, i.e. to 6.667 x 10-^ and 3.333 x 10"^. When the BGK 
model is employed for computations of drop oscillations, it is stable only to 
the point where the viscosity is lowered to u — 1.6667 x 10~^. Thus, the 3D 



X = 



{2n+l)'^{piPgPlPg)2 



(98) 



2 2 R^ [npg + {n + l)pi] (pipi) 2 + (pgpg) 2 
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Fig. 4. Configurations of an oscillating drop at different times in lattice units; 
pi/pg = 4:,fii/fig = 4,1^1 = Vg = 1.6667 X 10^2^ n = 0.1. 

MRT multiphase flow model allows the viscosity to be lowered (or enhances 
the Reynolds number) by a factor of 5 in this compared to the BGK 

model employing the same underlying physical model, based on the kinetic 
theory of dense fluids [13,14]. 

Figure 5 shows the interface locations of the oscillating drops as a function of 
time for the three shear viscosities noted above when k = 0.1. As expected, 
when the viscosity is reduced, it takes longer for the drops to reach the equilib- 
rium shape by viscous dissipation. The computed time periods of oscillations 
Tlbe are 1201, 1173 and 1159 when the viscosities are v = 1.6667 x 10~^, 
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Fig. 5. Interface location of an oscillating drop as a function of time for different 
kinematic viscosities, vf, pi/pg = 4:,iJ,i/ fig = 4, k = 0.1. 

6.667 X 10~^ and 3.333 x 10~^, respectively. The corresponding analytical time 
periods Tanai are 1258, 1202 and 1173, respectively. Thus the maximum rel- 
ative error is 4.80%. Let us now decrease the surface tension parameter k to 
0.08 from 0.10. Decreasing the surface tension parameter reduces the surface 
tension of the drop. As a result, it is expected that the drop would take longer 
to complete a period of oscillation. 

Figure 6 shows the interface locations for the three shear viscosities noted 
above when k — 0.08. The computed time period, for example, when i/ — 
1.6667 X 10~^ is 1412 and the corresponding analytical value is 1348. Both 
these values are higher than those at the same value of viscosity for the higher 
surface tension parameter, which is consistent with expectation. To confirm 
the trend with surface tension, we performed computations by increasing k 
to 0.12 from 0.10. Figure 7 shows the interface locations when k, — 0.12. 
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Fig. 6. Interface location of an oscillating drop as a function of time for different 
kinematic viscosities, vf, pi/pg = 4:,iJ,i/ fig = 4, k = 0.08. 

When u = 1.6667 x 10~^, the computed time period is 1143, which agrees 
with the analytical value within 4.60%. As this value is lower than with the 
surface tension parameter of 0.1, we conclude that the computations reproduce 
variations with surface tension which are consistent with expectation from the 
analytical solution. The 3D MRT model is able to reproduce the time period 
of oscillations within 5%. 

The 3D MRT model remains as stable as in the previous cases even for more 
complex dynamical problems such as when drops undergo shearing and other 
types of motions. In a recent study, this model was applied to simulate col- 
lision of a pair of drops under different conditions [37]. In particular, when 
employed to study off-center collisions it is able to reproduce experimentally 
observed [48] complex interfacial shape changes due to such dynamical events 
as shearing, rotation, stretching deformation and breakup, while maintaining 
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Fig. 7. Interface location of an oscillating drop as a function of time for different 
kinematic viscosities, vf, pi/pg = 4:,iJ,i/ fig = 4, k = 0.12. 

an improved stability compared to the BGK model. Indeed, in some cases re- 
ported in Ref. [37] as much as an order of magnitude improvement in stability 
was observed by the use of MRT model in heu of the BGK model. 

6 Summciry and Conclusions 

In this paper, we develop 3D MRT LB models for multiphase flows. They 
employ the LBE with a generalized colhsion term together with forcing terms 
representing interfacial physics. By considering 3D lattice velocity models and 
a second-order discretization of the forcing terms, the LBE written in terms 
of the distribution function is transformed into an equivalent system repre- 
sented by a set of conserved and non-conserved moments. The conserved or 
hydrodynamic moments include the fluid density and momentum, while the 
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non-conserved or kinetic moments include heat fluxes and viscous stresses. 
The models are developed such that when the collision term in the LBE is 
written in moment space, the moments relax to their equilibrium values at 
rates that can be adjusted independently. 

By applying the Chapman-Enskog multiscale analysis to the MRT models, 
we show that they correctly recover the 3D hydrodynamical equations for 
multiphase flows in the continuum limit. The dynamical relationships between 
various moments and the forcing terms, including surface tension forces, are 
systematically derived. Some of the relaxation parameters in the collision term 
are shown to be related to the shear and bulk kinematic viscosities of the fluid. 
The ability to independently adjust the relaxation parameters enhances the 
stabihty of the MRT models. The models are evaluated for accuracy by solving 
some multiphase test problems. It is shown that the 3D MRT models verify 
the Laplace- Young relation for static drops to within 8%. Computations with 
MRT models of drop oscillations show that shear viscosities can be lowered by 
a factor of 5 when compared to the BGK model. The computed time period 
of oscillations agrees with the analytical solution to within 5%. The MRT 
model, when applied to more complex multiphase problems such as binary 
drop collisions, is found to be as stable as that for these simple canonical 
problems. The MRT approach developed in this work can be extended to LBE 
multiphase flow models such as those developed in Refs. [16,17] to handle high 
density ratio problems. Moreover, it also provides a more natural framework 
to extend the LBE for more general situations such as viscoelastic or thermal 
effects in multiphase flows in future investigations. 
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A Appendix. Simplified Transformation Matrices for Unit Lattice 
Spacing 



For unit lattice spacing and time steps, i.e. c — 1, the transformation matrix, 
Eq. (23), for the D3Q15 model simplifies to 
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and that for the D3Q19 model to 
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